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SUMMARY  (U> 

A  simplified  method  of  modelling  2  dimensional  flow  of  propellant  doughs  in 
a  kneading  disc  of  a  twin  screw  mix-extruder  (SME)  has  been  developed.  The 
2  dimensional  model  is  a  first  step  in  the  development  of  a  full  3  dimensional 
model  of  flow  in  screws  and  kneading  blocks  in  SMEs.  The  model  is  based  on 
the  FAN  method  of  Tadmor  as  modified  by  Szydlowski  et  al,  and  it  aims  to  use 
less  computational  effort  tnan  alternative  models  based  on  finite  element 
analyses.  Apart  from  one  model  for  non-Newtonian  fluids,  which  has  several 
serious  flaws,  the  existing  models  have  assumed  that  the  viscosity  of  the  fluid 
was  Newtonian.  The  model  presented  here  givesa  better  method  of  calculating 
the  flow  of  non-Newtonian  fluids  than  the  previous  method,  and  it  produces 
more  realistic  flow  velocity  profiles.  • 


©  Commonwealth  of  Australia 


Author's  address: 

Ordnance  Systems  Division 
Weapons  Systems  Research  Laboratory 
PO  Box  1700,  Salisbury 
South  Australia 

Requests  to:  Chief,  Ordnance  Systems  Division 


UNCLASSIFIED 

/'• 


WSRL-TM-10/90 


TABLE  OF  CONTENTS 

Page 

1.  INTRODUCTION  1 

2.  THEORY  2 

2.1  Approximation#!.  Newtonian  viscosity  3 

2.2  Approximation  #2.  Averaged  power  law  viscosity:  velocities  calculated 

from  Newtonian  viscosity  9 

2.3  Approximation  #3,  Averaged  power  law  viscosity:  velocities  calculated 

from  power  law  velocity  profile  -1 

2.4  Approximation  #4.  Full  solution  5 

3.  RESULTS  AND  DISCUSSION  6 

4.  CONCLUSIONS  9 

5.  ACKNOWLEDGEMENTS  9 

REFERENCES  10 

TABLE  1.  FLOW  RATES  AND  AVERAGE  VISCOSITIES  IN  THE  CENTRE 

ELEMENTS  OF  THE  FLOW  CHANNEL  AND  MAXIMUM  PRESSURES  8 

APPENDIX  I:  DERIVATION  OF  FINITE  DIFFERENCE  EQUATIONS  13 

LIST  OF  FIGURES 

1.  Intermeshing  kneading  blocks  '2 

2.  Geometry  of  opened  out  flow  channel  for  modelling  17 

3.  Circulating  flow  in  the  channel  18 

4  Quantities  involved  in  flow  in  elements  1  and  1+1  18 


5  Flow  diagram  of  the  solution  process  for  the  equations  of  flow 


19 


WSRL-TM-10/90 


6.  Pressure,  velocity  and  channel  profiles  calculated  from  Approximation  #1 

7.  Pressure,  velocity  and  channel  profiles  calculated  from  Approximation  #2 

8.  Pressure,  velocity  and  channel  profiles  calculated  from  Approximation  #3 

9.  Pressure,  velocity  and  channel  profiles  calculated  from  Approximation  #4 

1 0.  Schematic  representation  of  the  combination  of  drag  and  pressure  flows 

1 1.  Pressure,  velocity  and  channel  profiles  calculated  from  Approximation  #1 
with  the  average  viscosity  calculated  from  the  mean  shear  rate 

12.  Pressure,  velocity  and  channel  profiles  calculated  from  Approximation  #4 
with  a  flight  tip  gap  of  0.5  mm 

13.  Pressure,  velocity  and  channel  profiles  calculated  from  Approximation  #4 
with  a  flight  tip  gap  of  1.00  mm 


ocesalon  For _ 

:i;is  GRA&I 
OTIC  TAB 

Llionnour.ced 

justification- 


T 


By - - - 

Distribution/ _ 

Availability  Codoa 
|Avall  and/or 
Special 


1 


WSRL-TM-10/90 


1.  INTRODUCTION 

The  study  of  the  flow  behaviour  of  propellant  materials  in  twin  screw  mix-extruders  (SMEs)  has  recently 
received  considerable  attention  in  the  USA.  Much  of  this  work  has  been  done  by  a  group  headed  bv 
L).  Kalyon  at  the  Stevens  Institute,  NJ,  and  also  by  R.  Armstrong  at  MIT.  Some  of  this  work  has  been 
reported  to  the  Continuous  Mixer  and  Extruder  Users  Group  Meetings  held  annually  in  the  US(ref.l,2),  as 
well  as  in  the  open  literature(ref.3,4).  The  work  was  aimed  at  producing  detailed  and  accurate 
descriptions  of  the  flow  in  SMEs  of  a  limited  range  of  materials,  usually  inert  composite  propellant 
analogues.  The  computer  models  used  in  the  analyses  were  fully  3  dimensional  and  used  complex 
constitutive  equations.  The  models  were  of  such  complexity  that  they  had  to  be  run  on  supercomputers  and 
graphics  processors.  These  facilities  are  not  available  for  flow  modelling  at  WSRL,  where  studies  aie 
being  undertaken  on  the  manufacture  of  nitrocellulose  propellants  by  SME. 

A  need  exists  for  the  development  of  modelling  methods  which  can  be  used  with  ordinary  mainframe 
computers,  or  preferably,  with  the  new  generation  desktop  computers.  These  simplified  models  would  not 
provide  the  same  degree  of  detail  as  the  complex  models,  but  they  could  be  used  as  an  aid  in  experimental 
design  by  providing  a  guide  to  the  likely  effect  of  changes  of  various  parameters  on  processing  behaviour 

An  approximate  2  dimensional  finite  element  analysis  of  flow  of  Newtonian  fluids  in  a  screw  channel  was 
developed  by  Denson  and  Hwang(ref.5),  and  a  similar  approach  to  modelling  non-Newtonian  fluids  was 
taken  by  Lai-Fook  et  al(ref.6).  However,  these  finite  element  approaches  are  too  complex  and  consume 
too  much  computer  time  to  be  useful  in  modelling  propellant  flow  in  SMEs.  A  simplified  model  of 
Newtonian  flow  in  a  kneading  disc  was  described  by  Werner  and  Eise(ref.7),  but  no  details  of  the  solution 
method  were  given. 

Simple  models  of  flow  in  SMEs  based  on  the  FAN  method  of  Tadmor(ref.8>  have  been  developed  by 
Szydlowski,  Brzoskowski  and  White(ref.9).  A  basic  approximation  of  these  models  is  the  reduction  in 
the  dimensionality  of  flow  from  3  to  2,  which  results  in  a  considerable  reduction  in  complexity  and 
computing  time  without  a  large  loss  of  accuracy  in  the  calculation  of  gross  flow  characteristics. 

A  more  drastic  approximation  from  the  viewpoint  of  propellant  processing  was  the  assumption  of  a 
Newtonian  viscosity.  The  assumption  of  a  Newtonian  viscosity  implies  that  the  viscosity  is  a  constant 
independent  of  the  shear  rate,  and  hence  it  is  independent  of  the  velocity  distribution.  This  assumption 
considerably  reduces  the  complexity  of  the  equations  describing  flow,  because  they  are  transformed  from 
non-linear  to  linear  equations.  One  important  consequence  of  linear  behaviour  is  that  solutions  of  the 
equations  remain  solutions  when  multiplied  by  constant  factors,  and  the  sum  of  two  or  more  solutions  is  a 
solution.  This  allows  a  wide  variety  of  situations  to  be  analysed  by  the  solution  of  only  a  few  equations. 

Unfortunately,  strong  non-Newtonian  behaviour  is  shown  by  most  polymeric  materials,  inc'uJing 
nitrocellulose  (NO  based  propellants(ref.lO),  and  the  assumption  of  a  Newtonian  viscosity  is  not 
acceptable.  Some  modelling  of  the  flow  of  non-Newtonian  fluids  has  been  done  by  Szydlowski  and 
White(ref.ll)  using  the  so-called  power  law  equation  for  viscosity.  This  equation  gives  a  reasonably 
accurate  description  of  viscosity  in  many  cases,  but  the  way  the  equation  was  applied  n  the  model  does 
not  appear  to  give  the  full  accuracy  possible  (this  matter  is  discussed  later). 
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The  aim  of  the  present  program  of  work  is  to  build  on  the  results  of  Szydlowski  and  White  to  produce  a 
more  realistic  model  which  is  applicable  to  the  flow  of  nitrocellulose  propellant  materials  in  SMEs.  This 
paper  describes  a  first  step  in  this  direction,  viz,  the  modelling  of  2  dimensional  flow  in  a  single  kneading 
disc.  There  are  two  reasons  for  starting  with  this  simplified  approach.  Firstly,  it  is  more  efficient  to 
develop  the  required  numerical  procedures  using  a  simplified  problem,  and  secondly  it  is  easier  to  explore 
various  more  accurate  approximation  procedures  which  will  be  used  in  more  complex  models. 


2.  THEORY 

A  typical  pair  of  kneading  blocks  made  up  of  a  staggered  series  of  specially  shaped  discs  for  a  co-rotating 
twin  screw  extruder  is  illustrated  in  figure  1.  In  an  actual  extruder  these  blocks  are  set-up  in  intermeshing 
pairs  along  parallel  shafts  inside  a  figure  8  shaped  barrel.  However,  for  the  purposes  of  modelling,  only 
the  behaviour  of  a  single  kneading  block  in  a  cylindrical  barrel  will  be  considered.  The  main  consequence 
of  this  assumption  is  that  the  effect  of  flow  in  the  intermeshing  region  between  the  pairs  of  blocks  is 
ignored. 

The  analysis  of  the  simplified  screw  geometry  is  similar  to  the  method  of  analysis  of  a  single  screw- 
ex  truderf  ref.  12).  The  main  assumption  in  the  analysis  is  that  (conceptually)  the  screw  channel  is  opened 
out  so  that  the  barrel  wall  becomes  flat,  and  the  channel  in  the  kneading  block  resembles  the  shape  in 
figure  2.  A  further  assumption  is  that  the  kneading  block  remains  stationary  and  the  barrel  wall  moves 
in  the  direction  opposite  to  the  screw  shaft  rotation. 

The  analysis  of  twin  screw  SMEs  also  involves  some  differences  from  the  analysis  of  single  screw- 
extruders.  The  flow  channels  in  single  screws  are  relatively  shallow  and  can  be  approximated  by  flat 
plates,  but  in  the  case  of  self-wiping  twin  screws,  the  channels  are  relatively  deep  and  have  a  complex 
shape,  and  hence  a  more  sophisticated  modelling  method  must  be  used. 

An  extra  simplifying  assumption  is  made  in  modelling  flow  in  screws  or  kneading  blocks,  viz,  that  the 
flow  is  essentially  2  dimensional.  Obviously  there  will  be  some  circulating  flow  in  the  channel,  and  near 
the  tips  there  will  be  flow  in  the  xz  direction,  see  figure  3.  Neglect  of  this  flow  may  reduce  the  accuracy 
of  the  model  slightly,  and  also  it  may  contribute  to  the  convergence  problems  at  small  flight  tip  gaps  (to 
be  described  later). 

For  the  modelling  of  flow  in  a  single  kneading  disc  considered  in  this  paper,  the  dimensionality  of  the 
flow  can  be  further  reduced  to  1.  It  is  assumed  that  motion  in  the  radial  direction  is  negligible  and  that 
the  circumferential  velocity  does  not  vary  rapidly  in  the  circumferential  direction.  These  assumptions 
are  summarised  as 


V]  =  V2  =  0  ,  and 


jjvj 

dx3 


=  0 


(1) 


The  general  equations  of  motion  for  a  viscous  fluid  can  be  reduced  for  the  geometry  of  figure  2  and  the 
assumptions  in  equation  (1)  to 
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For  common  fluids  the  shear  stress,  032,  can  be  written  in  terms  of  a  generalised  viscosity,  p,  and  the  shear 
rate,  3V3/3X2,  as 


<332 


(3) 


For  the  case  where  p  is  a  constant,  equation  (2)  becomes 


3P  _  32V3 

3x3  ^  3x2“ 


(4) 


This  equation  provides  the  basis  of  the  description  of  flow  of  the  material,  and  the  general  method  of 
solution  of  it  is  based  on  the  FAN  method  of  Tadmor(ref.8),  as  modified  by  Szydlowski  et  aKref.9).  The 
channel  is  divided  into  series  of  N  elements  in  the  X3  direction,  indexed  with  the  variable  I,  as 
illustrated  in  figure  4.  The  pressure  depends  only  on  X3,  and  it  is  assumed  to  be  constant  within  each 
element  and  vary  only  between  elements.  There  are  no  sources  or  sinks  for  flow  in  the  channel,  so  the  total 
flow  into  any  element  must  equal  the  flow  out  of  that  element.  In  particular  for  the  l'th  element  the  flow 
in,  Q(l-l),  must  equal  the  flow  out,  Q(I).  The  solution  is  achieved  by  deriving  equations  for  the 
magnitudes  of  the  flows  in  terms  of  the  pressures,  and  then  solving  the  equations  to  give  the  pressure 
distribution.  Once  the  pressure  distribution  has  been  determined,  it  can  be  used  to  calculate  the  velocity 
distribution. 

Details  of  the  methods  of  modelling  the  flow  are  given  below  in  order  of  improving  approximation  and 
increasing  expected  accuracy. 

2.1  Approximation#!.  Newtonian  viscosity 

The  simplest  equations  result  when  the  viscosity,  p,  is  taken  to  be  a  constant  which  is  independent  of 
shear  rate.  This  approximation  was  made  by  Szydlowski  et  al  in  reference  9.  In  this  case 
equation  (4)  is  linear,  and  can  be  integrated  directly  to  give  the  velocity  as 


V3<x2>  =  v3 . 2  -  sr  ■  r-  •  I(x2/h)-(x2/h)2] 

H  Z\i  dx3 


(5) 
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where  H  is  the  channel  depth  and  V3  is  the  velocity  of  the  wail.  Integrating  over  the  cross-section 
gives  the  volumetric  flow  rate: 


Q  = 


H  H3  3P 

3  '  2  '  12p  ■  ax3 


(6) 


Equation  (6)  can  be  applied  to  each  element  in  the  flow  channel  to  give  a  series  of  equations  for  the 
pressure  profile.  Details  of  the  equations  and  the  method  of  solution  are  given  in  Appendix  1.  Once 
the  pressure  distribution  has  been  determined,  it  can  be  substituted  into  equation  (5)  to  calculate  the 
velocity  distribution. 

2.2  Approximation  #2.  Averaged  power  law  viscosity:  velocities  calculated  from  Newtonian 
viscosity 

The  assumption  of  a  shear  rate  independent  viscosity  is  not  realistic  for  propellant  materials,  which 
are  known  to  show  considerable  shear  thinning  behaviour(ref.lO).  An  adequate  description  of  the 
viscosity  of  propellants  is  given  by  the  power  law. 


ft 


dv3  (n-1) 
3x2 


(7) 


where  k  is  the  consistency  index,  and  n  is  the  power  law  exponent  with  a  typical  value  of  0.4  for 
propellants.  Since  the  viscosity  varies  with  the  shear  rate,  which  in  turn  varies  with  depth  in  the 
channel,  the  viscosity  is  a  function  of  X2,  ie  p  =  gUj). 

Because  g  is  no  longer  a  constant,  equation  (4)  is  not  valid.  For  variable  values  of  g,  substituting 
equation  (3)  into  equation  (2)  gives 


3P  _  3g_  3v3  <^3 

3x3  3x2  3x2  3x22 


Substituting  equation  (7)  into  equation  (8)  gives,  after  some  manipulation 


3P  3^V3 

3x3  =  n 


32V3 

3x22 


(n-1) 


(9) 
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Hence  for  the  particular  case  of  a  power  law  dependence  of  viscosity,  the  Newtonian  equation  is 
modified  by  only  a  multiplicative  factor.  In  the  model  the  value  of  k  was  multiplied  by  n  and  the 
equations  solved  in  the  same  way  as  equation  (4).  It  appears  that  this  factor  was  not  taken  into 
account  by  Szydlowski  and  White(ref.9),  which  further  reduces  the  accuracy  of  their  results. 

The  simplest  way  to  take  account  of  the  variable  viscosity  is  to  calculate  an  average  viscosity  for 
each  element  of  the  channel,  see  figure  2,  and  substitute  this  value  in  the  equation  for  flow  ratetref.ti). 
The  average  viscosity  is  calculated  by  dividing  the  channel  depth  into  a  series  of  M  levels,  indexed 
with  the  variable  j,  and  calculating  the  viscosity  at  each  level,  then  averaging.  There  are  a  number 
of  types  of  average  which  could  be  used.  Szydlowski  and  White  in  reference  9  calculate  an  average 
shear  rate  and  then  substituted  it  into  equation  (7).  A  second  possibility  is  to  calculate  the  mean 
viscosity  from  the  values  in  each  layer,  but  this  method  was  found  to  cause  the  iterative  procedures  in 
the  model  to  diverge.  A  third  method,  which  was  found  to  be  the  most  accurate,  is  to  take  the 
reciprocal  of  the  mean  of  the  reciprocals  of  the  individual  viscosities,  ie 


f*av 


_ 1 _ 

average  (1  /pi) 


This  form  of  average  was  suggested  by  the  fact  that  the  viscosity  appears  in  the  denominator  of 
equation  (6). 

The  shear  rate  in  layer  J  is  determined  from  the  calculated  Newtonian  velocity  profile.  This 
approach  was  used  by  Szydlowski  and  White(ref.U). 

2.3  Approximation  #3.  Averaged  power  law  viscosity:  velocities  calculated  from  power  law 
velocity  profile 

In  reverse  pitch  kneading  blocks  and  screws  there  will  be  considerable  backflows  and  correspondingly 
large  velocity  gradients.  The  velocity  profile  calculated  from  the  above  procedure,  and  the 
corresponding  viscosities,  are  unlikely  to  be  related  to  the  true  values.  Hence  a  superior  method  from 
that  used  by  Szydlowski  and  White  is  required. 

As  a  first  step,  a  better  approximation  to  the  real  velocity  profile  can  be  calculated  from  equation  (9), 
using  the  already  calculated  viscosity  and  pressure  profiles  from  Approximation  42,  see  a  flow 
diagram  of  the  solution  process  given  in  figure  5.  The  method  of  solution  of  equation  (9)  is  given  in 
Appendix  1. 

2.4  Approximation  #4.  Full  solution 

The  pressures  and  volume  flow  rates  obtained  in  Approximation  #3  are  not  correct,  as  they  were 
calculated  using  averaged  viscosities.  The  full  solution  to  equation  (9)  must  be  obtained  by  an 
iterative  procedure  involving  adjusting  the  pressure  profile  to  balance  the  actual  inflow  and  outflow 
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of  each  element  with  the  volume  flow  rates  calculated  from  the  most  recently  calculated  velocity 
profile.  The  calculation  of  the  true  velocity  profile  must  also  use  the  individual  viscosities  from  each 
level,  and  the  process  is  illustrated  in  the  flow  diagram,  figure  5. 

The  computer  program  developed  in  this  study  contains  an  algorithm  for  varying  the  pressure  profile 
in  an  iterative  process  until  the  values  of  the  pressure,  velocity  and  viscosity  converge.  As  a  check 
that  the  solution  had  converged  to  the  correct  values,  the  pressures  were  also  calculated  directly  from 
equation  (9),  using  the  final  values  of  velocity  and  viscosity.  Details  of  the  calculation  are  given  in 
Appendix  1.  The  agreement  was  found  to  be  better  than  1%. 


3.  RESULTS  AND  DISCUSSION 

The  analysis  described  above  is  applied  to  modelling  the  flow  in  a  kneading  disc  in  the 
VVerner-Pfleiderer  ZSK  40  SME  at  WSRL.  The  barrel  diameter  is  40.3  mm,  and  the  distance  between 
screw  shafts  is  33.4  mm.  The  nominal  screw  element  diameter  is  40.0  mm,  which  gives  a  screw  tip-barrel 
clearance  of  0.15  mm.  Because  of  numerical  instability  problems  the  minimum  tip  clearance  modelled  was 
0.2  mm.  The  shape  of  the  disc,  and  hence  the  flow  channel,  was  calculated  by  the  method  of 
Booytref.13). 

The  power  law  parameters  for  viscosity  (see  equation  (7))  used  in  the  model  were  k  =  12500  Pa  and  n  =  0.4. 
These  values  are  typical  of  solvent  processed  double  base  doughs! ref.  10).  These  doughs  usually  show  a 
yield  stress  when  extruded  in  the  batch  process,  but  it  is  felt  that  the  structure  responsible  for  the  yield 
stress  would  be  broken  down  in  the  SME,  and  so  the  yield  stress  would  be  negligible. 

The  calculated  velocity  and  pressure  profiles  for  a  screw  speed  of  120  RPM  and  a  flight  bp  gap  of  0.2  mm, 
for  the  various  approximations,  are  given  in  figures  6  to  9.  The  bottom  part  of  the  figures  is  a  scaled 
representation  of  the  shape  of  the  flow  channel,  with  the  first  horizontal  line  representing  the  barrel 
wall,  The  middle  part  of  the  figures  shows  the  velocity  profiles  at  various  positions  along  the  length  of 
the  flow  channel.  The  dashed  lines  are  zero  velocity  for  the  appropriate  profile.  The  top  part  of  the 
figures  is  the  pressure  profile  in  the  channel. 

Only  the  leading  half  of  the  flow  channel  is  shown  in  the  figures,  as  the  profiles  are  symmetrical  about 
the  centre  of  the  channel.  The  pressure  profile  shows  that  the  pressure  built  up  by  the  flow  in  the  channel 
drops  linearly  across  the  flight  tip  to  the  maximum  negative  level,  and  then  builds  up  again  in  the 
channel.  In  the  trailing  half  of  the  channel  (not  shown)  the  pressure  reaches  the  maximum  positive 
pressure  just  before  the  flight  tip. 

The  calculated  pressure  profiles  are  qualitatively  similar  to  the  profiles  obtained  by  Szydlowski  and 
White(ref.9),  and  also  Werner  and  Eise(ref.7). 

The  velocity  profiles  are  complex,  and  show  the  effect  of  the  interaction  of  drag  flow  and  pressure  flow. 
For  pure  drag  flow  between  parallel  plates  the  velocity  profile  is  linear,  and  for  pure  pressure  flow  it  is 
parabolic.  The  addition  of  pressure  flows  to  drag  flows  produces  composite  profiles  due  to  the  presence  of 
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backflows,  see  figure  10.  It  can  be  seen  from  figures  6  to  9  that  the  velocity  profile  in  the  gap  is  convex  due 
to  the  pressure  gradient  forcing  the  material  forward  through  the  gap.  In  the  body  of  the  channel  the 
direction  of  the  pressure  gradient  is  reversed,  and  there  is  backflow  near  the  channel  root. 

The  Newtonian  profiles  are  given  in  figure  6.  The  maximum  pressure  is  30.8  MPa,  and  the  velocity 
profiles  in  the  channel  are  strongly  curved.  The  non-Newtonian  Approximation  #2,  see  figure  7,  shows  a 
much  smaller  value  of  maximum  pressure  of  1.17  MPa,  because  the  viscosity  is  greatly  reduced.  The 
velocity  profiles  are  unchanged  in  this  approximation.  This  case  corresponds  to  the  approximation  used 
by  Szydlowski  and  White(ref.U). 

In  Approximation  #3,  see  figure  8,  the  pressures  have  still  been  calculated  using  the  Newtonian  formula 
for  flow  rate,  but  the  velocity  profiles  have  been  calculated  using  a  variable  viscosity  profile,  and  it  can 
be  seen  that  the  velocity  profiles  are  considerably  flattened. 

The  full  solution  is  given  in  figure  9.  There  are  significant  differences  in  pressure  and  velocity  profile 
from  the  previous  approximations.  The  pressure  is  increased  by  about  50%  to  1 .65  MPa,  the  velocity 
profiles  are  further  flattened,  and  the  backflow  is  considerably  increased. 

The  calculated  flow  rates  in  the  flight  tip  gap  region  were  fairly  independent  of  the  approximation  used, 
because  the  flow  is  almost  totally  due  to  drag.  However,  in  the  centre  of  the  flow  channel,  the  calculated 
flow  rates  varied  considerably  depending  on  the  method  of  calculation,  see  Table  1.  It  is  readily 
apparent  that  the  total  flow  rate  from  the  central  elements  is  only  a  small  fraction  of  the  flow  occurring 
inside  the  central  elements,  and  the  total  flow  is  the  difference  between  large  forward  and  reverse  flows. 

For  the  Newtonian  case  the  flow  rates  were  calculated  directly  from  equation  (6),  and  also  by  integrating 
the  calculated  velocity  profile  using  Simpson's  Rule.  No  difference  could  be  found  between  the  two 
results.  The  change  in  calculated  flow  rates  between  Approximations  #1  and  #2  was  relatively  small 

However  there  was  a  significant  difference  between  the  calculated  flow  rates  using  Approximation  #3.  In 
this  approximation  the  pressure  profile  is  calculated  from  flow  rates  given  by  equation  (6)  by  substituting 
an  average  viscosity,  and  these  values  of  the  flow  rate  in  each  element  are  all  equal  to  each  other. 
However,  the  true  value  of  the  flow  rate  in  each  element  was  calculated  from  the  actual  velocity  profile, 
and  these  values  vary  between  2.79  x  10‘5  m3/s  in  the  flight  tip  gap  to  1.51  x  10"*  m3/s  in  the  centre  of  the 
channel.  This  indicates  that  Approximation  #3  does  not  give  a  very  accurate  estimate  of  the  pressures  or 
flow  rates,  and  it  is  not  usable  as  an  improved  approximation  over  Approximation  #2. 

It  was  possible  that  the  poor  results  from  Approximation  #3  were  due  to  the  method  of  averaging  the 
viscosity.  Szydlowski  and  White  used  an  average  viscosity  obtained  by  substituting  the  mean  shear  rate 
into  the  power  law  equation  (7).  Running  the  model  with  this  form  for  the  average  viscosity  gave  a  flow 
rate  for  the  centre  element  of  2.40  x  10“*  m3/s,  which  is  considerably  more  in  error  than  the  method  used 
in  this  work.  The  corresponding  pressure  and  velocity  profiles  are  given  in  figure  11,  and  comparing  with 
figures  8  and  9,  it  can  be  seen  that  the  velocity  profiles  are  further  from  the  true  profiles.  For 
Approximation  #2,  the  difference  in  flow  rates  calculated  from  the  different  averages  were  not 
significant,  but  the  pressure  calculated  from  viscosities  using  the  mean  shear  rate  was  about  7%  lower. 


WSRL-TM-10/90 


8 


Comparing  the  results  for  Approximations  #2  and  #4  in  Table  1,  it  can  be  seen  that  the  difference  in  total 
now  rate  is  only  of  the  order  of  5%.  However,  there  is  a  significant  difference  in  the  magnitudes  of  the 
forward  and  reverse  flows,  and  the  maximum  pressures  vary  by  about  50%.  The  velocity  profiles  in 
figures  6  and  9  show  a  large  difference  in  flow  patterns  and  hence  mixing  effectiveness.  These  results 
indicate  that  in  the  geometry  considered  in  this  work  the  extra  computational  effort  required  by  the  full 
solution  is  probably  justified. 

The  effect  of  varying  the  flight  tip  gap  can  be  seen  in  figures  9  and  12  to  13.  As  the  gap  widens,  the 
maximum  pressure  decreases  because  of  the  smaller  resistance  to  flow  over  the  flight  tip.  The  lower 
pressure  reduces  the  amount  of  flattening  of  the  velocity  profiles.  The  curvature  of  the  second  velocity 
profile  changes  from  convex  to  concave.  The  flow  rates  in  the  centre  elements  for  different  tip  gaps 
calculated  from  approximation  #4  are  given  in  Table  1.  The  total  flow  rate  increases  by  a  factor  of  4  for 
an  increase  in  tip  gap  of  a  factor  of  5.  Both  the  forward  flow  increases  and  the  back  flow  decreases  with 
increasing  tip  gap. 


TABLE  1.  FLOW  RATES  AND  AVERAGE  VISCOSITIES  IN  THE  CENTRE  ELEMENTS 
OF  THE  FLOW  CHANNEL  AND  MAXIMUM  PRESSURES 


Flow(m3/s  x  105) 

Approximation 

Total 

Forward 

Reverse 

Viscosity 

(Pa) 

Max  Pressure 

(MPa) 

#1 

2.62 

26.40 

23.80 

5000 

30.80 

#2 

2.79 

26.40 

23.80 

466 

1.17 

#3 

15.10 

23.20 

8.03 

625 

1.52 

#4 

2.90 

19.40 

16.50 

577 

1.65 

#3 

24.00 

27.30 

3.34 

511 

1.31 

(Average  |i  from 
average  shear  rate) 

Tip  Gap  0.5  mm 

#4 

8.22 

21.70 

13.50 

611 

0.96 

Tip  Gap  1.0  mm 

#4 

11.80 

26.40 

8.36 

670 

0.56 

Study  of  the  velocity  profiles  calculated  using  Approximation  #4  indicates  where  the  maximum  stresses 
are  developed,  and  where  viscous  heat  generation  would  be  greatest.  In  the  flight  tip  gap  the  site  of 
maximum  shear  and  heating  is  near  the  tip  and  away  from  the  wall,  in  the  channel  the  shear  is  highest 
near  the  barrel,  and  this  has  a  shorter  path  to  the  barrel  wall,  and  so  should  aid  heat  transfer  from  the 
propellant  material  to  the  circulating  coolant  in  the  barrel. 
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It  can  be  seen  that  plug  flow  occurs  in  the  lower  part  of  the  channel,  and  mixing  is  confined  to  the  upper 
layers  in  the  channel.  This  type  of  flow  has  been  seen  in  flow  visualisation  studies  of  kneading 
blocks(ref.3).  Forward  and  reverse  kneading  blocks  with  a  staggering  angle  of  30°  showed  unmixed  areas 
in  the  centre  of  the  flow.  Flow  in  these  blocks  would  resemble  the  situation  modelled  here,  where  the 
staggering  angle  is  effectively  0°. 


4.  CONCLUSIONS 

A  simplified  method  of  modelling  the  2  dimensional  flow  of  propellant  materials  in  a  kneading  disc  of  a 
SME  has  been  developed.  The  model  uses  more  accurate  approximation  methods  than  used  previously  for 
calculating  the  non-Newtonian  flow  shown  by  propellant  materials.  The  model  also  includes  a 
multiplicative  factor  which  was  omitted  from  previous  work.  The  method  should  be  applicable  to 
modelling  the  3  dimensional  flow  in  kneading  blocks  and  conveying  screws. 
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APPENDIX  I 

DERIVATION  OF  FINITE  DIFFERENCE  EQUATIONS 
1.1  Pressure  profile  in  the  Newtonian  approximation 

The  flow  rate  Q  in  the  X2  direction  of  a  Newtonian  fluid  with  viscosity  p  between  2  parallel  plates  a 
distance  H  apart  and  moving  with  a  relative  velocity  Vj  is  given  by 


Q  =  V3 


Hi  iL 

1 2r  3x3 


(6) 


where  3P/3x2  is  an  applied  pressure  gradient. 

In  the  flow  channel  defined  by  figure  2  there  are  no  sources  or  sinks  for  flow,  so  the  total  flow  into  any 
element  must  equal  the  flow  out  of  the  element  In  particular  for  the  l'th  element  the  flow  in,  Qil-I). 
must  equal  the  flow  out,  Q(I),  see  figure  4.  The  flows  are  evaluated  at  the  edges  of  the  elements,  so 
the  variables  associated  with  the  flow  must  be  evaluated  at  the  same  points  These  variables 
include  H([),  V3U). 

The  differential  equation  (6)  is  solved  by  converting  it  to  a  finite  difference  equation.  Separate 
equations  are  set  up  for  flows  into  and  out  of  each  element,  and  these  are  equated.  Hence 


Q(l-l)  =  Q(l> 


(1.1) 


and 


H(l-l)  H(l-l)3  3PU-1) 

2  12p  3xs 


HU)  H(l)3  dm 

2  12q  3x2 


Converting  this  equation  to  a  difference  equation  gives 


(H(l)  i-  H(l-ll)  V3  UHU)  s  HU-l  ))/2)3  (P(t)  •  PU-D)  7 

2  2  12q(l>  '  DX(l)  -  DX(l-l)  ‘ 


(H(ls-l)  ■»  H(D)  V3  «H(M)  +  HU))/2)3  (P(ls-l)  -  PCD) 

2  '2  12gU+1)  DXU+U+  DXU)  '  2 
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Since  the  number  ot  equations  is  relatively  small,  they  can  be  solved  by  an  efficient  direct  substitution 
technique! ref.  14).  The  boundary  conditions  are  determined  by  the  cyclic  nature  of  the  geometry.  The 
value  of  the  pressure  in  the  O'th  element  is  the  same  as  for  the  N  +  1st  element,  and  when  these 
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conditions  are  fed  into  the  equations,  a  solution  is  possible.  However  the  absolute  magnitude  of  the 
pressure  is  not  determined,  and  so  a  constant  is  added  to  the  calculated  pressures  to  give  zero  pressure 
in  the  center  of  the  flight  tip. 

1.2  Solution  of  the  differential  equation  for  velocity 

The  velocity  profile  can  be  calculated  from  the  differential  equation  describing  flow, 

^  =  n.k.(^2|(n'1).^2  (9) 

8x3  t  3x2  J  3x2“ 


where  3p/3x3  is  known.  This  equation  is  a  second  order  differential  equation,  and  can  be  readily 
solved  by  finite  difference  methods. 

The  difference  equation  is  derived  in  the  following  way.  The  first  derivatives  are  given  by 


3v3<X2) 

3x2 


V3(|-l>-  v3(|) 
dxs 


for  X2  >  xs(J) 


and 


3V3U2) 

3x2 


V3U)  -  V3(j-1) 

dx2 


for  x2  <  xs<J ) 


where  dx2  is  the  step  in  the  X2  direction  between  X2  -()•’■  I  )/H  and  xi  =  I  /  H.  The  step  size  in  uniform  in 
this  case.  The  second  derivative  of  V3  is  the  derivative  of  the  first  derivatives. 


32V3(X2) 

3x22 


V3U  +  I)  -  V3U )  V3CI)  -  V3U-I) 

dxs  dxs 

— - fc - - - 

dx2 


=  1V3(]^1)  -2v,(|)  *  V3()-l)|/dx22 


(1  41 


For  the  J'th  level  in  the  1'th  element  equation  (9)  becomes 


VVSRL-TM-10/90 


16 


v3(J + 1 )  -  2v3(J)  +  v3(J- 1 ) 


=  dx2*  . 


|P(I  +  1)  -  P(l)l 
[DXll+1)  +  DX(I)]/2 


(1.5) 


Those  equations  have  the  tridiagonal  form  of  equations  (1.6),  and  they  are  solved  in  the  same  way. 
I J  Calculation  of  the  pressure  profile  from  the  velocity  profile 
Equation  (9)  is  used  to  calculate  the  pressure  profile,  ie. 


3P 

3x3 


(9) 


The  corresponding  difference  equation  is  a  rewritten  form  of  equation  (1.6) 


P(U1)-P(I) 


[v3(H)-2v3(J)  -  v3(H)) 
dx32 


p.(l,j>  .  IDXd+l)  *  DX(l)|/2 


Since  the  values  of  velocity  v3(l)  and  viscosity  pd)  are  known  from  the  previous  iteration,  the  values 
of  the  pressure  increments  between  elements,  P(I+1)  -  Pd),  can  be  calculated.  The  pressure  profile 
P(x3)  is  obtained  by  successively  adding  P(l+1)  -  P(I)  to  the  previously  calculated  pressure  P(l)  to  give 
the  ordinate,  and  adding  (DX(I+1>  +  DX(l))/2  to  the  previously  calculated  x3  to  give  the  abscissa. 
The  starting  point  is 


(P(l)),x3(0»  =  (0,0). 
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Figure  2.  Geometry  of  opened  out  flow  channel  for  modelling 


Figure  3.  Circulating  flow  in  the  channel 


Figure  4.  Quantities  involved  in  flow  in  elements  I  and  1+1 
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APPROXIMATION  #4 


Figure  5.  Row  diagram  of  the  solution  process  for  the  equations  of  flow 
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Figure  6.  Pressure,  velocity  and  channel  profiles  calculated  from  Approximation  #1 


Figure  7.  Pressure,  velocity  and  channel  profiles  calculated  from  Approximation  #2 
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Figure  8.  Pressure,  velocity  and  channel  profiles  calculated  from  Approximation  #3 


Figure  9.  Pressure,  velocity  and  channel  profiles  calculated  from  Approximation  #4 
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Drag  * 
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Drag  ♦ 
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Figure  U). 


Schematic  representation 


of  the  combination  of  drag  and  pressure  flows 
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Max  Prommuram  965  IMra)  Parameter* 

K-  1 2500  (Pa) 

N-  .  4 


VELOCITY  PROFILES  Hex  Velocity*  -251  <in/e> 


Figure  12.  Pressure,  velocity  and  channel  profiles  calculated  from  Approximation  #4 
with  a  flight  tip  gap  of  0.5  mm 


Power  Law 

Max  Preeeure-  -.55 1 000000000000 1  (MPa)  Paramo ter* 

K-  12500  (Pa) 

N«  .  4 


VELOCITY  PROF  I LE5  M0«  V. 1 oc l ty-  . 251  <«/.) 


CHANNEL  DEPTH  PROFILE 


Figure  13.  Pressure,  velocity  and  channel  profiles  calculated  from  Approximation  #4 
with  a  flight  tip  gap  of  1.00  mm 
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